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ABSTRACT 

Most problems in gravitational lensing require numerical solutions. The most frequent types of problems are (1) 
finding multiple images of a single source and classifying the images according to their properties like magnifica- 
tion or distortion; (2) propagating light rays through large cosmological simulations; and (3) reconstructing mass 
distributions from their tidal field. This lecture describes methods for solving such problems. Emphasis is put on 
using adaptive-grid methods for finding images, issues of spatial resolution and reliability of statistics for weak 
lensing by large-scale structures, and methodical questions related to shear-inversion techniques. 



1. Introduction 



Only for very special lens models can numerical methods be 
avoided in gravitational lensing studies. There are three essential 
^ reasons for that. One is the non-linearity of gravitational lensing, 
. the fact that image and source positions are related to one an- 
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other in a non-linear fashion. This gives rise to the well-known 
phenomena of mutiple imaging, strong image distortions, and 
so forth. The second reason is that lenses exist which are them- 
selves best described by numerical models. Galaxy clusters are 
one example, lensing by large-scale structures is another. Al- 
though it is true that many aspects of gravitational lensing by 
^3 large-scale structures can be derived analytically, detailed sim- 
Qnulations require numerical techniques. The third reason is that 
the interpretation of gravitational lensing effects or events often 
require the application of sophisticated algorithms to ever grow- 
ing amounts of data. One example is the reconstruction of the 
projected mass density distribution of a galaxy cluster from the 
observed image distortions due to gravitational shear. 
. ^ Needless to say, there are many more aspects of numerical 
K> methods related to gravitational lensing than I can cover in this 
^ review. An outstanding example are the highly elaborate meth- 
od ods that have been developed over recent years for determining 
" image shapes of faint background galaxies on CCD frames, and 
for extracting the gravitational shear signal from them. This is a 
whole branch of data analysis on its own. Here, I can only deal 
with numerical methods for relating mass distributions to their 
gravitational lensing effects. 

Consequently, the outline of this lecture is as follows: First, 
I shall discuss methods for studying individual lenses, i.e. their 
imaging properties, their critical curves and caustics. In par- 
ticular, the use of adaptive grids and techniques for searching 
and characterising images will be discussed. Second, I shall de- 
scribe how extended lenses can be treated numerically using the 
multiple-lens plane theory. This will lead to the basic equations 
for tracing light rays through (simulated) cosmological volumes. 
A large fraction of the discussion will be devoted to issues of res- 
olution and noise, and to spurious effects in simulated lensing 
statistics. Finally, third, I shall describe inversion techniques, 
i.e. methods for reconstructing the projected mass distribution of 
lenses whose distortion has been measured. The classic Kaiser- 
Squires method will be described, and also maximum-likelihood 
techniques and maximum-entropy methods. 

General lensing theory and the theory of weak lensing are 
covered by Koenraad Kuijken's and Peter Schneider's lectures 
in this volume. Basic references on lensing include the textbook 
by Schneider et al. (1992) and the lecture by Narayan & Bartel- 



mann (1999), reviews of weak lensing are Mellier (1999) and 
Bartelmann & Schneider (2001). 



2. Individual Lenses 

2.7. Assumptions 

A brief reminder of the basic assumptions underlying the the- 
ory of individual lenses may be in order. There are three main 
assumptions. First, the Newtonian gravitational potential of the 
lens be small, |<I>| <C c^. Second, velocities in the gravitational 
lens system, both of constituents within the lenses and of the 
lenses with respect to the rest frame of the microwave back- 
ground, be small v <C c. Third, the extent of the lenses along 
the line-of- sight be small compared to the other distances in the 
system, which are usually cosmological and thus comparable to 
the Hubble radius, c/Ho = 3h~^ Gpc, with Hq being the Hubble 
constant and /z = //q / 1 00 km s~ ^ Mpc~ ^ . 

It is worth noting how well these assumptions are satisfied in 
ordinary lensing situations. Consider a galaxy cluster with mass 
M = 10^^ h~^MQ. Assuming spherical symmetry, the Newto- 
nian potential at a distance R= lh~^ Mpc from its centre is 
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evidently much smaller than the speed of light squared. A 
typical length scale for the radius of a galaxy cluster is 1 — 
1.5 Mpc, which is several hundred times smaller than typical 
distances in a cluster-lensing system. Finally, peculiar velocities 
of galaxy clusters with respect to the Hubble flow are of order 
several hundred kms~^, and typical velocities of galaxies within 
galaxy clusters reach of order 10^ kms~^ but both velocities are 
way below the speed of light. The above assumptions hold even 
better for lensing by galaxies, of course. 

We can thus safely assume the above conditions to be satis- 
fied. It is then possible to project the lensing mass distribution 
onto a plane perpendicular to the line-of-sight, the lens plane, 
and describe it by its surface mass density Z. Sources are as- 
sumed to be located on a corresponding plane, the source plane. 
A typical lens system is sketched in Fig. [I] 

The three distances /)d,s,ds shown in Fig. [I] and explained in 
its caption are generally not additive because of space-time cur- 
vature, thus Z)s 7^ Z)d + ^ds in contrast to flat space-time. 



2.3. The Lensing Potential 
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Fig. 1. — Schematic view of a gravitational lens system. The lens is 
projected onto the lens plane perpendicular to the line-of- sight, sources 
are located on the parallel source plane. There are three distances re- 
quired to describe the geometry of the system, i.e. the distances i)d,s,ds 
between the observer and the lens, the observer and the source, and 
between the lens and the source, respectively. Due to space-time curva- 
ture, these distances are generally not additive. 



2. 2. Coordinates and Notation 

Let us now introduce physical coordinates ^ and f| on the lens 
and source planes, respectively. Alternatively, it is often conve- 
nient to introduce angular coordinates and P, which are obvi- 
ously related to % and x\ through 



Dimensional coordinates are of course not suitable for numerical 
calculations, which can only handle numbers. We thus have to 
introduce a length scale ^o, or alternatively an angular scale 0o, 
in the lens plane. This length scale is so far arbitrary. It implies 
a length or angular scale 
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in the source plane. Dimension-less coordinates are then defined 
by 
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in the lens and source planes, respectively. The numerical code 
will have to deal with the dimension-less vectors x and y. It helps 
numerical accuracy greatly if these numbers are of order unity. 
Thus, the first challenge in setting up a lensing simulation is to 
choose an appropriate length- or angular scale or 0o, which 
should both be adapted to the physical problem at hand, and to 
the requirement that numerical codes work most accurately if the 
numbers they are dealing with are neither too large nor too small, 
compared to machine accuracy. Choosing unappropriate length 
scales can, for instance, render image searches unsuccessful. 



It will be convenient for the following discussion to introduce 
the lensing potential \|/ as the basic physical quantity for lens- 
ing studies. It is the scaled, projected Newtonian gravitational 
potential of the lens. 



(5) 



The so-called reduced (i.e. appropriately scaled) deflection an- 
gle is the gradient of the potential. 



a{x) = V3e\i/(x) , 



(6) 



and the lensing convergence (i.e. the appropriately scaled 
surface-mass density) is 
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(7) 



Finally, the gravitational tidal field is described by the two- 
component shear. 



yi (x) = - - \i/,22) = ^ (ai,i - a2,2) , 72 (^) = ¥,12 = oci,2 , 

(8) 

where the convention was used that fij is the derivative of the i- 

th component of / with respect to the coordinate xj. It is impor- 
tant to note that the fact that all lensing quantities can be derived 
from the scalar lensing potential establishes relations between 
all of them. This will be exploited several times later. 

Note that the lensing quantities must be rescaled in case the 
coordinate scale ^0 is changed. Suppose is introduced instead 
of ^0- Since the physical surface-mass density of the lens must 
remain the same at any given physical location, the reduced de- 
flection angle must transform as 



a(/) = ^a(jc), 
so 



(2) and convergence and shear transform as 



[k(x'),yK^')] = (|) [k(x),y,-(x)] . 



(9) 
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2.4. Imaging 

Suppose now we were given some description of the lensing po- 
tential or of the deflection angle a{x). This description 
could be an analytical formula, or it could be in form of an ar- 
ray, i.e. a set of numbers given at grid points (xi^xj). We wish to 
know how the given lens images its background. 

We introduce a coordinate grid Xfj on the lens plane subject 
to the condition that it be sufficiently well resolved. This means 
that the smallest features in the lens must be covered by at least 
a few grid points. Since we are given the deflection angle as 
a function of position, we can compute a deflection-angle grid, 
OLij = OL{xij). The mapped grid on the source plane is then simply 
ytj = Xij — OLij. This mapped grid will appear as a distorted image 
of the regulargrid in the lens plane, as the example in the left 
panel of Fig. |2| shows. 

The mapping process must now be reversed in order to obtain 
an image created by the lens. For doing so, the source plane 
is first covered with a regular grid, f^j. Next, we loop over all 
grid points Xij in the lens plane and find its mapped source point 



Fig. 2. — Left panel: A regular grid in the lens plane (blue dots) is 
mapped onto the source plane (red dots) using a numerical description 
of a deflection- angle field. Distortions are clearly visible. Right panel: 
For each point in the lens plane, those points of a regular grid in the 
source plane (blue) are searched which surround its mapped point in 
the source plane (red). 



jij in the source plane, and search for the nearest neightbours 
surrounding yji in the source plane. This is illustrated in the right 
panel of Fig.|2| The surface brightness of the source, known at 
the positions j^^, can then be interpolated to jij and the result 
assigned to the image point That way, the surface brightness 
at all points in the lens plane can be determined, and thus the 
lensed image be constructed. Fig. O shows an example. 

The left panel of the figure shows a simulated CMB tempera- 
ture fluctuation field of 10' x 10' angular size. The temperature 
increases from white to red. In essence, the temperature fluctua- 
tion corresponds to a fairly smooth gradient across the field. The 
right panel shows the gravitational lensing signature imprinted 
on the CMB at such angular scales by a galaxy cluster. The tem- 
perature visible at an angular position on the sky, r'(0), is re- 
lated to the intrinsic temperature T through r'(0) = r [0 — a(0)] . 
Thus, the light deflection by the cluster causes the visible tem- 
perature distribution to be rearranged, yielding a highly specific 
pattern (Seljak & Zaldarriaga 2000). 

2.5. Critical Curves and Caustics 

As mentioned in the introduction, the deflection- angle field con- 
tains full information on the lensing mass distribution. All other 
quantities like convergence and shear, but also image magnifica- 
tions, follow from the deflection angle via differentiation. It is 
thus a common task to compute numerical derivatives. 

Suppose a function fix) is tabulated on a grid, so that we 
are given the values fij at the grid points Xij. The derivative of 
fix) at a particular point xqo in the first coordinate direction is 
approximated by 

^ =:^(/io-/-io) + o(/z^) (11) 

where h is the separation of the grid points in the chosen direc- 
tion; cf. the left panel of Fig.|3 This centred difference has the 
advantage compared to the more straightforward one-sided dif- 
ferences /lo — /oo or /oo — /-lo of being second-order in the grid 
separation h. There are higher-order differencing schemes using 
function values at more than two grid points, but the second- 
order scheme is usually sufficient. No lensing quantity should 
vary strongly between two adjacent grid points because other- 
wise the resolution of the grid would be grossly insufficient. 

We will typically need derivatives of the deflection angle a. 
Since a is itself the gradient of a scalar potential, its derivatives 



Fig. 4. — Left panel: Second-order numerical differentiation using cen- 
tred differences. Right panel: A simple method for finding points in the 
lens plane next to a critical curve uses sign changes of the Jacobian de- 
terminant between the point considered and its four nearest neighbours. 



must satisfy ai^2 = V,i2 = V,2i = 0C2,i- It is thus usually prefer- 
able to check that this relation is satisfied within numerical ac- 
curacy, and to use (ai^2 + 0C2,i)/2 instead of either ai^2 or a2,i 
alone. 

Critical curves in the lens plane are defined by the condition 
that the Jacobian determinant of the lens mapping vanish there, 
det (x) =0. The elements of the Jacobian matrix are j^^j = 
hij — aij, thus the Jacobian determinant is 

Z) = detJ^ = (l-ai,i)(l-a2,2)-a?,2 • (12) 

It can be computed once the (numerical) derivatives of the both 
deflection-angle components have been determined. 

One method of identifying grid points in the lens plane next 
to the critical curve proceeds as follows. Let S = sign(Z)), and 
consider one particular grid point xqq in the lens plane. The 
point is next to the critical curve if, and only if, the sign of the 
Jacobian determinant changes between it and one or more of its 
nearest neighbours. Hence, if the condition 

^00(^-10 + ^10 + ^0-1+%) < 4 (13) 

is satisfied, the gridpoint xqq is next to a critical curve (cf . the 
right panel of Fig. |4ll. Of course, xqo is not itself on the criti- 
cal curve, but to the positional accuracy determined by the grid 
resolution, the position of the critical curve can be constrained 
that way. Points on the source plane next to the caustic curve are 
then easily found via the lens equation, yaj = xaj — dL{xcij), 
where the xaj are the grid points in the lens plane next to criti- 
cal curves. 

As an example, consider a lens model for a spiral galaxy, con- 
sisting of a spherical halo and a flat disk seen almost edge-on 
(Bartelmann & Loeb 1998). The deflection-angle field of such 
a lens can be given analytically (cf. Keeton & Kochanek 1998). 
Convergence and total shear (yi +72)^^^ as determined by nu- 
merical differentiation are shown together with the modulus of 
the deflection angle in Fig. 13 

The critical curves and caustics of that lens model as deter- 
mined with the method described above are shown in Fig.|6| 

2.6. Adaptive Source Grids 

One of the most prominent goals of gravitational lensing studies 
with individual strong lenses is to determine the imaging statis- 
tics of a given lens model, for example the abundance of highly 
magnified events, the occurrence of multiple imaging with the 
images satisfying certain conditions, and the like. This is done in 
principle by distributing many sources across the source plane, 
imaging them as described before, and determining the image 




Fig. 3. — Left panel: A simulated CMB temperature fluctuation field of 10' x 10' size. Right panel: The same field, lensed by a galaxy cluster, 
which imprints a characteristic pattern on the temperature fluctuations. 
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Fig. 6. — Critical curves (left) and caustics (right) of the spiral-galaxy 
lens model illustrated in Fig.|5] 



properties. However, such events are rare. If one were to cover 
the entire source plane with a regular grid of sources, this grid 
would have to have a very high resolution for rare events to be 
reliably found. In turn, most of the sources probed would pro- 
duce images failing the criteria imposed, so by far the largest 
fraction of the CPU time used would be wasted. 

This situation calls for adaptive grids. We know in advance 
that any strongly lensed image will occur near a critical curve, 
or any strongly lensed source near a caustic. It is those sources 
that we need to treat in detail, while those far from caustic curves 
are usually only required to normalise the statistics properly. 

One approach for defining an adaptive grid, and there may be 
others more suitable for a particular lensing situation, proceeds 
as follows. Again, we assume that we know the deflection an- 
gle of the lens, either because it was provided numerically or 
because it is described by a known analytic formula. Then, we 
saw in the preceding subsection how grid points can easily be 
identified which are close to a critical curve in the lens plane, or 
a caustic curve in the source plane. 

In order to save computational time, the source plane is first 



covered with a coarse grid. This grid should obviously be fine 
enough for the caustics to be properly resolved; for instance, 
it must not be so coarse that the typically two types of caustic 
curve, the radial and the tangential one, are closer than a few 
times the grid separation. 

Next, those points on that coarse grid are identified and saved 
which are next to a caustic curve. This can, for instance, be 
done by masking, i.e. by attaching a logical variable to all grid 
points and setting it to either true or false depending on whether 
it is or is not next to a critical curve. One can then cover the 
source plane with a grid whose resolution is doubled in both 
dimensions, and keep only those points which are identical with, 
or surrounded by, points of the coarse grid which were masked 
in the preceding step. This procedure can be repeated as often 
as desired, i.e. until the finest grid level reaches the ultimately 
required resolution. Note that it is not the grids and their masks 
which need to be saved, but only the coordinates of those grid 
points which are either part of the coarse initial grid, or whose 
logical mask values are true. That way, lists of source positions 
can be constructed which are to be probed later for the images 
they give rise to. 

Naturally, this can only be a basic recipe which needs to be 
adapted to the situation at hand. For instance, the condition that 
grid points need to be next to a caustic can be replaced by the 
condition that the absolute value of the Jacobian determinant be 
less than a given threshold which can be lowered at each step 
of grid refinement. Such a criterion would naturally increase 
the grid resolution near such grid positions where sources are 
certain to be highly magnified. 

Of course, if statistics is the ultimate goal, one has to take 
into account that sources near caustic curves were positioned 
such as to have an unfair advantage over sources far from caus- 
tics. Since we have chosen to double the grid resolution at each 
refinement step, each source on a refined grid represents only a 
quarter of the area on the source plane represented by a source 
on the next coarser grid. Assigning a statistical weight of unity 
to the sources on the finest grid, the weight must quadruple for 



Fig. 5. — A lens model for an almost edge-on spiral galaxy: shown are the modulus of the deflection angle (left), the convergence (centre), and the 
absolute value of the shear (right). 



each coarser level. If the grid was refined N times, the weight 
of sources on the coarsest grid is thus Wi = 2^^. Each source is 
assigned a statistical weight Wi in that way, and counts Wi times 
in the final statistical evaluation. 

The left panel in Fig. 7 shows the source locations chosen for 
evaluating image statistics of the spiral lens model illustrated in 

Figii 

2. 7. Finding Images 

The principle of finding the images of a given source is simple: 
Given a source at position %, find those grid points Xij on the 
lens plane which are mapped sufficiently close to %, i.e. whose 
mapped points yij are within a specified distance from % . 

The problem with this approach is that a square- shaped or 
rectangular grid cell from the image plane is mapped onto a dis- 
torted figure in the source plane. In most cases, this figure will 
be a parallelogram, but in rare cases, opposing corners of the 
original rectangle may even be interchanged on the source plane. 
How can it then be decided whether a given point in the source 
plane is inside or outside the mapped grid cell, or in other words, 
whether the image of the given source falls within that particular 
grid cell on the lens plane? 

The solution is to split each grid cell in the lens plane into 
two triangles, because a mapped triangle always remains a tri- 
angle, which always has a well-defined interiour (cf. Schneider 
et al. 1992). 

Consider Fig. [HI The three grid points marked on the lens 
plane in the left panel of the figure are mapped to the distorted 
triangle shown on the right panel, which contains the source po- 
sition. Call <ii,2,3 the three vectors from the mapped triangle's 
corners towards the source position. It can be shown that the 
source is inside the mapped triangle if the three vector products 



X J2 , X ^3 , J2 X (14) 

are all positive, with the vector product in two dimensions being 
defined as 

dxb = a\b2 — a2b\ . (15) 

One straightforward way to verify this condition is to con- 
vince one's self that the source point is inside the triangle if 
all vectors di point within the angles spanned by the adjacent 
sides of the triangle, and that this condition translates to Eq. (IT4ll 
above. 




Fig. 8. — Illustration of the technique for finding images described in 
the text. Grid cells in the lens plane are split into triangles (left panel), 
which have a well-defined interior after being mapped back onto the 
source plane (right panel). This would not necessarily be the case for 
rectangular grid cells. A source is contained by a triangle if all mixed 
cross products di x dj for the shown vectors di are positive. 



This algorithm for finding images works well as long as the 
separation between images is larger than the size of the grid cells 
in the lens plane. Very close images can be contained within the 
same grid cell, in which case the algorithm would find only one. 
Of course, this potential problem can be remedied by increasing 
the grid resolution, but then a very large number of grid cells 
would have to be checked in vain for containing an image. 

Again, a viable solution uses adaptive grids. One can start 
with a coarse grid on the lens plane. Searching for images on 
that coarse grid will almost certainly not yield all images of a 
multiply imaged source, but those missed will be closer than the 
grid separation to those found. Then, those grid cells containing 
images can individually be covered with a highly resolved grid, 
and the image search repeated on those sub-grids. Hence, the 
first step represents a coarse scan of the lens plane for grid cells 
containing at least one image, and the second step scans only 
those regions on the lens plane in detail where images are sure 
to be found. If needed, further sub-grids can be similarly nested. 

Of course, even though this procedure is highly adaptive and 
efficient, it always has a remaining resolution limit, and images 
closer than that will not be resolved. It is then important to adapt 
the resolution of the finest sub-grid to the situation at hand, for 
instance such that remaining unresolved images would neither 
be resolved by observations. The right panel in Fig. illus- 
trates the result of an adaptive image search for all sources at the 
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Fig. 7. — ^L^/r panel: source positions placed on a multiply refined grid in the source plane. Caustic points are black. Obviously, the source plane 
is best sampled near caustics. Right panel: Number of images found for sources placed at the positions shown in the left panel. 



positions shown in the figure's left panel. Colours denote im- 
age numbers: Black means one image, blue three, and red five, 
while green shows source positions for which an even number of 
images has been found, in contradiction to the necessarily odd 
image number produced by non- singular lenses. Such events are 
rare, but they do occur because of the finite resolution limit of 
the algorithm applied. 

Figure|9lgives an example for possible results of that adaptive 
technique for finding images. Colour-coded is the total magni- 
fication of point sources in the source plane behind the almost 
edge-on spiral lensing galaxy introduced above. The increas- 
ing spatial resolution towards the caustic curves is evident. The 
panel inserted into the figure shows caustics (blue) and critical 
curves (red) of the lens, the source position as a blue dot just 
inside the right-hand "naked" cusp, and the three images as red 
hexagons whose size logarithmically encodes the image magni- 
fication. 

2.8. Asymmetric Lenses 

So far, we have used a model for a spiral galaxy as an exam- 
ple for a complex lens whose properties need to be determined 
numerically. Despite its complexity, the model is still highly 
symmetric; and what is more, its deflection angle is given as an 
analytic formula. Sources were so far assumed to be point-like. 

Let us now increase the level of complexity and use a numer- 
ically simulated galaxy cluster to gravitationally lens extended 
sources. Again, we assume the deflection angle to be given and 
postpone the question as to how it can be determined from an 
A^-body simulation. 

All techniques described above for computing convergence 
and shear from the deflection angle, for finding critical curves 
and caustics, for placing sources on an adaptive grid, and for 
finding images within grid cells split into triangles remain valid 
unchanged. Figure [lO| shows an example. 

The modulus of the cluster's deflection angle is shown as the 
colour plot in the left panel. The right panel shows a section 




Fig. 9. — The colour encodes the total magnification of a point source 
lensed by an almost edge-on spiral galaxy; blue means a magnification 
near unity, yellow means very high magnification. The adaptive reso- 
lution of grid cells on the source plane is clearly visible. The size of 
the grid cells decreases substantially towards regions of high magni- 
fication. The inserted panel shows caustics and critical curves of the 
same lens (blue and red, respectively), a source position close to the 
right-hand "naked" cusp, and the three images as red hexagons, whose 
size logarithmically encodes their magnification. 



of the source plane with the dots marking source positions, and 
their colour illustrating the image number. Black, blue and red 
means one, three, or five images, respectively. The caustic struc- 
tures can clearly be identified as the boundaries between black 
and blue and between blue and red, respectively. 

2. 9. Imaging Extended Sources 

Extended sources can be described in a variety of ways. What 
follows is a simple description for elliptical sources, but alterna- 
tive source models can easily be constructed along similar lines. 
We assume that source positions % have already been found, 
preferentially on an adaptive grid as described before. Also, we 
need to be sure that the grid resolution in the source plane is suf- 





Fig. 10. — The colour plot in the left panel shows the modulus of the deflection- angle field in the lens plane of a numerically simulated galaxy 
cluster. The right panel shows how sources are adaptively placed on the source plane (dots), and how many images these sources have; one (black), 
three (blue) or five (red). The boundaries between the colours mark the caustic structure. 



ficiently high as to resolve the smallest sources to be considered. 

Elliptical sources are described by three more parameters, 
viz. their size, their ellipticity, and their position angle (|). Let us 
describe the ellipticity hy e = b / a, with a and b being the semi- 
major and semi-minor axes of the ellipse, respectively. Finally, 
we introduce an effective radius r by demanding that a circle 
of radius r have the same area as the ellipse, hence r = ^/ab. 
By rotating by an angle an ellipse centred on the coordinate 
origin whose axes are aligned with the coordinate axes, it can 
straightforwardly be shown that a grid point yij is enclosed by 
the ellipse if the condition 

+ 25ji5j2 sin (|) cos (|) ^- — <0%) 

is satisfied, where 83; = y^j If the grid point Xij, whose 

image in the source plane is yij, satisfies Eq. (IT6l) . the image 
point Xij is part of the source, and the image can be constructed 
by assigning the source's surface brightness at yij to the image 
point Xij. By mapping the entire lens plane onto the source plane 
and checking Eq. (IT6l) for each individual imaged grid point yij, 
all image points belonging to the given source can be identified. 

It is often desired for statistical purposes to automatically 
characterise a large number of images. An example is the de- 
termination of cross sections for the formation of large gravita- 
tional arcs by a numerically simulated galaxy cluster, for which 
a large number of sources need to be imaged and the image prop- 
erties automatically quantified to search for the rare "giant" arcs. 
Most of the methods described here have been introduced and 
used extensively e.g. by Bartelmann & Weiss (1994), Bartel- 
mann et al. (1995, 1998), Meneghetti et al. (2000, 2001); see 
also the contribution by Massimo Meneghetti to this volume. 

A source may have multiple images, thus the point sets in 
the lens plane found by imaging extended sources need not be 



connected. The first step is therefore to group the image points 
into images. This can be done with a variant of the classical 
friends-of-friends algorithm: Pick one arbitrary point out of any 
given set of image points and search for another image point 
which is at most V^h grid units away from the first point; h is 
the grid size in the lens plane. If there is such a point, it is called 
a "friend" and grouped into the same image as the first point. 
Now take the "friend" and repeat until no further "friends" can 
be found and the image is complete. If more image points are 
left on the lens plane, pick one of those and repeat the process 
until all image points have been grouped. If the image is large 
enough, and the grid resolution on the lens plane is high enough 
for the image to consist of many points, the image magnification 
is simply the ratio between the numbers of pixels covered by an 
image and the number of pixels covered by the source. 

Once all image points belonging to a single image have been 
identified, it is often useful to determine the boundary points 
of that image, e.g. by identifying those points inside an image 
which have a neighbour outside the image. By suitably order- 
ing the boundary points, a boundary line can be found whose 
length can be measured and used in further steps of the auto- 
matic image classification. Next, the curvature of the image 
can be found by first identifying the image point which of the 
source centre, then search for the boundary point most distant 
from the so-defined image centre, and finally searching for the 
boundary point most distant from the first boundary point. These 
three points uniquely define a circle whose radius can be used as 
an approximation for the arc radius. And so on, you get the 
drift: Once image points are grouped into individual images and 
boundary curves have been determined, images can be classified 
by adapting elementary geometrical figures to them. 

2. 1 0. Deflection Angles of Asymmetric Lenses 

So far we have assumed to be given the deflection angle either 
as an analytic expression or as two two-dimensional arrays of 



numbers giving its two components as a function of position in 
the lens plane. We now need to describe methods for obtaining 
the deflection angle of a numerically simulated lens. 

The first issue to be discussed is the spatial resolution. Since 
the simulated lens is composed of discrete particles which rep- 
resent a smooth mass distribution in reality, the deflection angle 
must not be computed by simply summing up the deflection an- 
gles of the individual particles: The result would be a collection 
of microlenses rather than a single macrolens, having many spu- 
rious and undesired imaging features. 

Rather, the collection of particles has to be projected onto a 
lens plane, on which it needs to be smoothed in some way. We 
will return later to the issue of how particles should be sorted 
into grid cells. An important point to be addressed before is 
how large the grid cells should be chosen. They should be 
small enough for important features of the lens to remain identi- 
fiable; they should be large enough for the surface density to lose 
the "graininess" due to its being composed of individual parti- 
cles, and they should be large enough so that Poisson errors are 
smaller than a certain threshold. If the number of particles per 
grid cell is nh^, its Poisson fluctuation is Vnh^, thus the discrete- 
ness of the particles gives rise to fluctuations in the surface-mass 
density. Demanding that the relative fluctuations of the density 
should be smaller than 8 <C 1, the cell size h has to be chosen 
such as to satisfy {nh^)~^^^ < 8. It is impossible to give a gen- 
eral rule applicable to the majority of lensing situations, but it 
is clear that resolution, smoothing and particle noise have to be 
carefully balanced by choosing the grid cell size appropriately. 

Assigning particle masses to grid points in order to obtain a 
smooth density distribution is an art of its own (cf. Hockney 
& Eastwood 1988). In principle, the particle mass could sim- 
ply be attributed to the single grid point next to its position. 
This "nearest grid point" (NGP) method is appropriate for par- 
ticles near the centre of a cell, but particles near cell boundaries 
should be attributed to the cell and its neighbour(s) in order to 
avoid boundary effects like density discontinuities. Numerous 
schemes for interpolating particles across cells have been pro- 
posed. They are generally of the form 



(17) 



where Q is the quantity to be interpolated onto a point x, e.g. the 
particle mass, the sum extends over all particles sufficiently 
close to the point of interest x, and W{x — Xi) is a smoothing 
or interpolation kernel depending on the separation vector be- 
tween the particle position Xi and x. The kernel is decomposed 
into three factors directions. 



W(5x) = w{5xi)w{5x2)w{5x3) , 



(18) 



one for each dimension, the i-th of which depends only on the 
/-component of the separation vector. Interpolation methods can 
now be classified according to the kernel factors w(5x) and their 
width. 

The "cloud-in-cell" (CIC) scheme uses the kernel factors 



wcic(8x) = |j '^-^'/'^ Otherwise" ' ^^^^ 



\5x\/h for \5x\ < h 



which implies that the particle is distributed over the four nearest 
grid points. A more elaborate scheme is the "triangular shaped 
cloud" (TSC) method, which uses the kernel factors 

3/4-dx^/h^ for \dx\<h/2 

WTSc(^) = { {3/2-\5x\/hf/2 for h/2 < \5x\ < 3h/2 . 
otherwise 

(20) 
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Fig. 11. — The "cloud-in-cell" (left panel) and "triangular shaped 
cloud" (right panel) interpolation schemes are illustrated here. The 
(projected) particle position is marked red, the grid points to which the 
particle mass is assigned are marked blue and green. The CIC and TSC 
schemes assign the particle mass to the eight and 27 nearest neighbours, 
respectively (in three dimensions). 



The CIC and TSC interpolation schemes are illustrated for two 
dimensions in Fig. [H] For all schemes, the kernel has to be 
normalised such that all particle mass fractions add up to unity. 

Suppose now we have obtained the surface mass density on a 
grid Kij = k(x;j), then the deflection angle can most straightfor- 
wardly be determined by direct summation as 



a/ 



^ij ^kl 
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(21) 



Depending on the number of grid cells, the direct summation can 
be prohibitively slow. In many circumstances of astrophysical 
interest, fast-Fourier techniques can then be applied. In order to 
see how this works, note that the deflection angle can be written 
as a convolution of the convergence k(x) with a kernel 



K{x) 



1 X 



(22) 



This allows the Fourier convolution theorem to be applied, 
which holds that the Fourier transform of a convolution is the 
product of the Fourier transforms of the functions to be con- 
volved, hence 

d.{k) = k{k)k{k) . (23) 

The Fourier transform of the kernel K can be determined and 
tabulated once. Using fast-Fourier techniques to determine the 
Fourier transform of the convergence k(^) requires the conver- 
gence to be periodic on the lens plane. In many cases, this 
can be safely assumed or arranged. Often, lens planes are con- 
structed from large-scale A^-body simulations which have peri- 
odic boundary conditions by design, or the lens is an isolated 
object like a galaxy cluster, which can be surrounded by a suf- 
ficiently large field for the convergence to drop near zero every- 
where around the edges of the field. Fast-Fourier methods speed 
up the computation of the deflection angle considerably. 

If necessary, derivatives of the deflection angle field can also 
be determined in Fourier space. Once the convergence has been 
Fourier transformed, one can employ the two-dimensional Pois- 
son equation to compute the Fourier transform of the lensing 
potential, 

2 . 



(24) 



from which the Fourier transforms of the deflection angle and 
the shear components can easily be determined. 
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V, y2 = -^i^2V- (25) 



Relations like those and the exploitation of fast-Fourier methods 
are particularly relevant for simulating gravitational lensing by 
large-scale structures. 

3. Lensing by Large-Scale Structures 

3. 1 . Resolution Issues 

Obviously, the thin-lens approximation that we have been using 
so far breaks down if one wishes to study gravitational lensing 
by large-scale structures. The solution then is to cover the com- 
plete cosmic volume whose lensing effects one wants to sim- 
ulate with simulation boxes stacked along the line-of- sight, to 
project suitable slices on individual lens planes, and to use mul- 
tiple lens-plane theory for describing light propagation. 

The multiplicity of lens planes, and the general weakness of 
lensing by large-scale structures, make questions of angular and 
mass resolution particularly relevant for cosmic lensing. For in- 
stance, lens planes close to the observer are typically poorly re- 
solved because even small grid cells span a large solid angle near 
the observer, and making grid cells smaller is not generally an 
acceptable solution because then the number of particles per grid 
cell becomes small, and the shot noise possibly unacceptably 
large. However, lens planes near the observer are less efficient 
than lens planes approximately half-way to the source because 
the lensing efficiency function is zero at the observer and source 
redshifts and peaks in between. Yet, structures grow over time, 
thus the lensing efficiency function is skewed towards lower red- 
shifts because structures are geometrically less efficient lenses, 
but their density contrast keeps growing. By a related argument, 
sources at very high redshifts do not require the entire cosmo- 
logical volume between them and the observer to be filled with 
lens planes because lens planes at very high redshift are geomet- 
rically inefficient and have a low density contrast. The left panel 
of Fig. El shows two examples for the lensing efficiency func- 
tion times the linear growth factor, which is the relevant quantity 
combining structure growth with geometrical efficiency. 

Similarly, the effective angular resolution of the simulation is 
dominated by the angular resolution of those lens planes near 
the peak in the combined efficiency function, i.e. the product of 
geometrical efficiency and linear growth factor. 

The shot noise caused by the discretisation of mass into par- 
ticles is particularly important for studies of weak lensing by 
large-scale structures. Even in absence of density inhomo- 
geneities, shot noise leads to density fluctuations. They need 
to be sufficiently smaller than the signal, i.e. the convergence 
fluctuations which cause weak lensing. 

In essence, this requirement also imposes a resolution limit. 
Suppose we wish to quantify the weak-lensing signal within a 
solid angle 5^1. The volume spanned by 5^1 within redshifts z 
and z + dz is 



da 



prop 



dz 



dz , 



(26) 



where D{z) and Dprop{z) are the angular diameter and proper 
distances to redshift z. In absence of density inhomogeneities, 
this volume element contains dN{z) particles, with 



dN{z)=dV{z)^, 



Mr, 



(27) 



where p(z) is the mean matter density at redshift z, and mp is the 
mass of an A/^-body particle in the simulation. The contribution 
to the lensing convergence by these particles has to be weighted 
by the effective lensing distance, Z)eff(z,Zs), and by numerical 



factors. Poisson fluctuations in the particle number thus cause 
convergence fluctuations whose variance is 



5^Koc r dzDif^{z,Zs)dN{z) . 
Jo 



(28) 



These fluctuations need to be compared with, and smaller than, 
the convergence fluctuations due to large-scale structure, which 
are typically of order (k^)^/^ ^5% for sources near redshift 
unity and angular scales of order V. According to Eqs. (l26l 
through (l28l . the rms shot noise (5^k)^/^ scales like 5^1^/^, thus 
the requirement that the signal-to-noise ratio 
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1/2 



(29) 



exceed a specified threshold translates into a lower limit to the 
solid angle 5Q. which can reasonably be resolved by the sim- 
ulation. The smallness of the rms cosmic convergence Knns = 
^^2^1/2 ijYipiies that many particles need to be enclosed by the 
"cone" spanned b y 3Q . for the simulation to be reliable. The 
right panel of Fig.TT2l shows an example. The rms cosmic con- 
vergence in per cent and the noise-to-signal ratio are plotted as 
functions of angular scale. The noise level was adapted to an N- 
body simulation with particle mass mp = 6.8 x 10^^ Mq. The 
curves show that the noise-to-signal ratio drops below unity for 
sources at redshift Zs = 1 only if the angular resolution is lowered 
to > 5', while an angular resolution of > 0.8' can be achieved for 
Zs = 1000 (i.e. for weak lensing of the CMB; Pfrommer 2002). 

3.2. Multiple Lens-Plane Theory 

Weak lensing by large-scale structures requires the cosmic vol- 
ume to be split into multiple lens planes rather than a single one 
(for general reference on multiple lens-plane theory, see Schnei- 
der et al. 1992). The lens plane closest to the observer is the 
image plane which represents the observer's sky. A light ray 

piercing the image plane at a physical coordinate is mutiply 
deflected on N lens planes and finally reaches the source plane 
at the physical coordinate 



Ti(^i) = -^^i + XA-sa(^0 



(30) 



where the Df and Di^ are the angular diameter distances from the 
observer to the /-the lens plane, and from the i-th lens plane to 

the source, respectively. The light ray passes the i-th plane at 4/, 

where it is deflected by a(^;). Similarly, the are determined 

Iy(li) = ^|i + IA7S(I), (31) 



i=l 



where Dij is the angular diameter distance from the i-th to the 
7-th lens plane. 

Introducing angular coordinates 0; = t,i/Di yields 



e,(eo = ei + '£^a(e,), 



(32) 



where we have introduced the reduced deflection angle a = 
{Dis/Ds)dL. We now define the matrices 
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(33) 
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Fig. 12. — Left panel: The product of lensing efficiency function times the Hnear growth factor for density perturbations is shown for two different 
source redshifts, Zs = 3 and Zs = 1000, respectively, the latter being relevant for gravitational lensing of the cosmic microwave background. The 
growth factor skews the geometrical lensing efficiency towards lower redshifts. At redshift 5, the combined efficiency function drops to 10% of its 
peak value for Zs = 1000, implying that by far not the complete redshift range up to Zs needs to be covered with lens planes. Right panel: The solid 
curves show the rms cosmic convergence for sources at three different redshifts in per cent, the dashed curves the noise-to- signal ratio obtained in 
an N-hody experiment with particle mass mp = 6.8 x \0^^h-^MQ. Both types of curve are plotted as functions of angular scale. 



Clearly, is the Jacobian matrix of the lens mapping between 
the /-th lens plane and the image plane, thus j^^v is the Jacobian 
matrix of the mapping between the source and image planes. 
The goal is thus to determine J^^v in order to obtain convergence, 
shear, and magnification for a light ray starting out into direction 
9i . The ray-tracing equation (l32l implies the recursion relation 



4. Inversion Techniques 

Let us conclude with a brief discussion of inversion techniques. 
They are typically less demanding numerically, but the methods 
which have been developed for this purpose are interesting in 
their own right. 



/ - > — ^ Ui^i , 



(34) 



starting with J^i = /, the identity matrix. In summary, the 
deflection-angle fields OC/ on the N lens planes can be used to 
construct the matrices Ui according to Eq. (l33t . then Eq. (l34l) 
can be used to determine the lensing experienced by a light ray 
starting out into any direction on the image plane. The left panel 
in Fig. [21 shows the total convergence experienced by sources 
at Zs = 5 on a lens plane with a side length of 4.25°, obtained 
from an A/^-body simulation (Pfrommer 2002). 

The right panel in Fig. [El shows numerically determined 
power spectra for the effective convergence as functions of wave 
number /, which is the Fourier conjugate variable to the angular 
scale. The lines in this figure show the theoretically expected 
power spectra. The agreement between the numerical and theo- 
retical results is very good over a limited range of wave numbers. 
Once the wave numbers increase beyond the limit set by the an- 
gular resolution, the simulated convergence fields lack power 
and the numerical results fall below the theoretical ones. This 
happens at lower / for smaller source redshifts, because a fixed 
angular scale, and thus wave number /, corresponds to smaller 
physical scales at lower distances. On the low-/ end, i.e. for 
large structures, the errors on the numerically determined power 
spectra increase because the number of independent modes in 
the simulated convergence field decreases as the modes increase. 
This example should suffice to demonstrate that numerical sim- 
ulations of gravitational lensing by large-scale structures should 
be carefully designed to match their final purpose. 



4. 1 . Shear Deconvolution 

We have seen before in Eqs. (l24ll and (l25ll that convergence and 
shear are related because they are both linear combinations of 
second derivatives of the scalar lensing potential \|/. In Fourier 
space, the relations are algebraic and can easily be combined to 
eliminate the Fourier transform \j/ of the potential. Transforming 
back into configuration space, the convergence turns out to be a 
convolution of the shear with a well-known kernel. 
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with 
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This is the classic Kaiser & Squires (1993) shear inversion equa- 
tion. Its limitations have been discussed in detail and removed to 
satisfaction by modifying e.g. the kernel components 2),^; they 
are not of interest for the discussion here (cf. Peter Schneider's 
lecture in this volume). 

A suitable practical approximation of using measured 
galaxy ellipticities 8,- (/ = 1 , 2) is 



K(6): 



1 ^ 

— E[®i^i.' + ®2e2,i] 



(37) 



where n is the number density of lensed galaxies on the sky. In 
practice, however, it turns out that an approximation like dTTl 
would have infinite noise because of the random sampling of 
the shear components by N galaxy ellipticities 8^. This can 



Fig. 13. — Left panel: Effective convergence on a field of 4.25° side length for sources at redshift Zs = 5, obtained using multiple lens plane theory 
on an N-bo&y simulation. Right panel: Effective-convergence power spectra measured with the same set of simulations (crosses), compared with 
theoretical expectations (lines), for different source redshifts. The numerical results follow the theoretical curves very well within an intermediate 
range of wave numbers /. At larger /, i.e. for small structures, the resolution limit of the simulation is reached and the power spectra fall rather 
steeply. At the low-/ end, the noise increases because the number of modes in the simulation decreases as the modes increase. The numerical 
power spectra for low-redshift sources fall below the theroretical expectation at lower / than for high-redshift sources because a given angle, and 
thus wave number /, corresponds to a larger physical scale at smaller distances from the observer (from Pfrommer 2002). 



be remedied by introducing a smoothed kernel (D ' instead of (D 
e.g. .2 2 

1 + ^ )exp(-^ 



0)' = 
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(38) 



where 0s is the angular smoothing scale (Seitz & Schneider 
1995). The noise convariance matrix between the convergence 
values at two different grid points 0; and dj is then 



K(eOK(0,-: 



exp 



202 



(39) 



where is the scatter of the intrinsic galaxy ellipticities (van 
Waerbeke 2000). This expression demonstrates that smoothing 
introduces correlations on the convergence map on the angular 
scale 0s, but the variance of k can become very high if 0s is 
chosen too small. A careful balance between the local variance 
and non-local correlations is necessary in order to arrive at a 
convergence map with the required properties. 

4. 2. Maximum-Likelihood Lens Inversion 

An entirely different approach to lens inversion uses the 
maximum-likelihood technique (Bartelmann et al. 1996). Each 
lensed background galaxy / provides a measurement of two el- 
lipticity components (£i,n£2,i) and its angular size. Comparing 
the size of a galaxy behind a galaxy cluster to the average size of 
unlensed galaxies of the same surface brightness, an estimate 
of the inverse magnification of the lensed galaxy can be derived. 
Thus N galaxies provide a 3A/^-dimensional data vector 



d = (ei,i,£2,bn,...,£i,iv,e2,iv,riv) . 



(40) 



The goal of the lens inversion is then to find a two-dimensional 
array \\fjk of lensing potential values such that the ellipticities 
and inverse magnifications caused by that potential at the posi- 
tions 0; of the real galaxies optimally reproduce the measured 
ellipticities and inverse magnifications. In other words, the po- 
tential values \\fjk have to be determined such as to minimise the 

mean- square difference between the data vector d and the model 
data vector d [\\fjk {xi ) ] , 



where the errors a, can be estimated from the data themselves. 
The minimisation of with respect to the potential values 
y\fjk can be done with any minimisation algorithm like, e.g. the 
downhill simplex method. For large fields, the number of po- 
tential values can become very large. In that case, conjugate- 
gradient methods are preferred, which make use of the fact that 
the derivatives of with respect to the \\fjk are known analyti- 
cally. Such methods can speed up the minimisation sufficiently 
to render it feasible even for large potential arrays (cf. Press et 
al. 1992). 



4.3. Maximum-Entropy Methods 

The minimisation of is a special case of the maximum- 
likelihood technique for assumed Gaussian deviations of the 
measured data around the model values. Improvements of the 
maximum-likelihood technique can be derived starting from 



Bayes' theorem, 
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which states that the probability P{"^\d) of finding the potential 

Xj/ given the data d is proportional to the probability P{d\\]f) of 
obtaining the data given the potential, times the probability /^(\|/) 

for finding the potential. The denominator P{d) is called the 
evidence and simply normalises Eq. (l42ll . /^(\|/) is called the 
prior, quantifying any a priori information one has or assumes 
on the potential \[/, P{d\\]f) is called the likelihood, and P(yf\d) is 
the posterior probability. The goal is now to maximise the latter, 
which is equivalent to maximising the product P{d\\]f)P{y]f) of 
likelihood and prior. If we have or can assume Gaussian noise 
and a diagonal noise correlation matrix, the likelihood reduces 
toP(J|v)=exp(-xV2). 

It can now be shown that in absence of any further informa- 
tion, the best, i.e. least prejudiced, prior is the maximum-entropy 
prior, 

P(yf) oc exp [aS(yf, m)] , (43) 



with the cross entropy 



3N 
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-\|//ln 
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(44) 



where m is a model vector for the potential which can encode 
expectations on the potential, or simply be chosen to be uniform 
for all /. The potential array is then determined by maximising 
exp(— %^/2 + a^"), or equivalently by minimising 
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instead of the simple in Eq. (l4Tl . The parameter a can be 
included into the minimisation. Bayesian theory implies that a 
good approximation to the optimal choice for a is determined 
such that F ^ 3N/2 at the potential minimum The error co- 
variance matrix for the potential \\f is given by the inverse curva- 
ture matrix of F, 



((\|/-\j/)(\|/-\j/)T); 



d^F 



(46) 



Maximum-entropy methods have been suggested and used for 
regularising shear-inversion techniques such that their spatial 
resolution is adapted to the strength of the lensing signal (Bridle 
et al. 1998; Seitz et al. 1998). 



5. Concluding Remarks 

Many numerical methods have been used for gravitational lens- 
ing studies which I was not able to cover during the limited time 
of the lecture. Among them are the hierarchical tree-code meth- 
ods introduced into microlensing by Wambsganss et al. (1990) 
and the methods for constraining cluster mass distributions from 
multiple arc systems (e.g. Kneib et al. 1993; see also Jean-Paul 
Kneib 's presentation in this volume). Despite this unavoidable 
incompleteness, I hope to have given a flavour of how numeri- 
cal methods can be used for lensing, and what the main problem 
areas are. 



